/**
* @file    ReMap.cpp
* @brief   Construct the "Move" sequence that remaps PropArray
*          Also change GridInfo, Gross2Net, Net2Gross
* @date    2021/08
* @version
* Copyright (c)
*/

#include <iostream>
#include <map>

#include "BinInput.h"
#include "SparseVector.h"

using namespace std;

const int MAX_STRING_LEN = 256;
const int REG_INFO_SIZE = 10;
const int iDead = -1;
const double NearZero = 1E-10;
const double InfSmall = 1E-20;

struct int_doub {
    int_doub() {};
    int_doub(const int _i, const double _db) {
        i = _i; db = _db;
    };

    int i;
    double db;
};

static int ia_to_str(const int* ia, const int MAX_SIZE, char* c) {
    int i = 0;
    for (i = 0; i < MAX_SIZE && ia[i] != 0; ++i) {
        c[i] = (char)ia[i];
    }
    int str_len = min(i, MAX_SIZE - 1);
    c[str_len] = '\0';

    return str_len;
};

static int ia_compare(const int* a1, const int* a2, const int len) {
    int different = 0;
    for (int i = 0; i < len; ++i) {
        if (a1[i] != a2[i]) ++different;
    }
    return different;
};

template<class _TYPE>
static bool Key_to_Ind(const std::map<_TYPE, int>& MAP1, const _TYPE& KEY, int& IND)
{
    typename std::map<_TYPE, int>::const_iterator itFind = MAP1.find(KEY);
    if (itFind != MAP1.end())
    {
        IND = itFind->second;
        return true;
    }
    else 
    {
        return false;
    }
};

int BinReader::ReMap(const double* GridInfo_in, const int* SubregInfo_in, const int* Host_map_in, const int* Gross2Net_in, const int* Net2Gross_in)
{
    Moves.clear();

    /* ��ȡӳ��ǰ������ṹ(Pre) */
    const UGINT64* Pre_GridInfo = origin_GridInfo;

    PropArray a_Pre_Subreg_info;
    PropArray a_Pre_Host_map;
    PropArray a_Pre_Porv;
    GetStaticData("SUBRINFO", a_Pre_Subreg_info);
    GetStaticData("HOSTMAP", a_Pre_Host_map);
    GetStaticData("PORV", a_Pre_Porv);

    /* ����쳣��� */
    // ���Pre_SubregInfo�ĳ���
    if (!a_Pre_Subreg_info.is_empty() && (a_Pre_Subreg_info.length < 2 || a_Pre_Subreg_info.length != a_Pre_Subreg_info[0] * a_Pre_Subreg_info[1] + 2)) {
        return -1;
    }
    // ���Pre_Host_map�ĳ���
    if (!a_Pre_Host_map.is_empty() && a_Pre_Host_map.length != (int)Pre_GridInfo[2]) {
        return -1; // ���Ȳ�������������
    }
    if (a_Pre_Host_map.is_empty() && (Pre_GridInfo[6] > 0 || Pre_GridInfo[7] > 0)) {
        return -1; // �м�������Ƕ��ʽ�ѷ�, ȴû��Host_map_in
    }
    // ���DIMENSƥ��
    if (Pre_GridInfo[0] < 3) {
        // CPG check NX, NY, NZ
        if ((int)GridInfo_in[0] == 3 || (int)GridInfo_in[3] != (int)Pre_GridInfo[3] || (int)GridInfo_in[4] != (int)Pre_GridInfo[4] || (int)GridInfo_in[5] != (int)Pre_GridInfo[5])
            return -2;
    }
    else if (Pre_GridInfo[0] == 3) {
        // GPG check NG, NLy
        if ((int)GridInfo_in[0] < 3 || (int)GridInfo_in[14] != (int)Pre_GridInfo[14] || (int)GridInfo_in[15] != (int)Pre_GridInfo[15])
            return -2;
    }
    // POROLVLƥ��
    if ((int)GridInfo_in[1] != (int)Pre_GridInfo[1]) {
        return -2;
    }
    // UNITƥ��
    //if ((int)GridInfo_in[19] != (int)Pre_GridInfo[19]) {
    //    return -2;
    //}
    /* ��ȡ����ṹ */
    int* Pre_Subreg_info = a_Pre_Subreg_info.is_empty() ? nullptr : a_Pre_Subreg_info.idata;
    int* Pre_Host_map = a_Pre_Host_map.is_empty() ? nullptr : a_Pre_Host_map.idata;

    // Gross sizes
    int PoroLvl = (int)Pre_GridInfo[1];
    int Pre_Gross_size = (int)Pre_GridInfo[2];
    int Pre_Gross_fin_size = (int)Pre_GridInfo[6];
    int Pre_Gross_emdf_size = (int)Pre_GridInfo[7];
    int Pre_Total_NLv0 = (Pre_Gross_size - Pre_Gross_emdf_size) / PoroLvl;

    int Gross_size = (int)GridInfo_in[2];
    int Gross_fin_size = (int)GridInfo_in[6];
    int Gross_emdf_size = (int)GridInfo_in[7];
    int Total_NLv0 = (Gross_size - Gross_emdf_size) / PoroLvl;

    int Base_NLv0 = (Gross_size - Gross_fin_size - Gross_emdf_size) / PoroLvl;

    // Net sizes
    int Pre_Net_size = (int)Pre_GridInfo[8];
    int Pre_FLvlBaseNum = (int)Pre_GridInfo[9];
    int Pre_FLvlFinNum = (int)Pre_GridInfo[10];
    int Pre_MLvlBaseNum = (int)Pre_GridInfo[11];
    int Pre_MLvlFinNum = (int)Pre_GridInfo[12];
    int Pre_emdf_size = (int)Pre_GridInfo[13];

    int Net_size = (int)GridInfo_in[8];
    int FLvlBaseNum = (int)GridInfo_in[9];
    int FLvlFinNum = (int)GridInfo_in[10];
    int MLvlBaseNum = (int)GridInfo_in[11];
    int MLvlFinNum = (int)GridInfo_in[12];
    int Emdf_size = (int)GridInfo_in[13];

    // Sub-region numbers
    int Pre_FinG_Num = 0;
    int Pre_EMDF_Num = 0;
    int record_size0 = 0;

    if (Pre_Subreg_info != nullptr) {
        int Pre_Num_Subreg = Pre_Subreg_info[0];
        record_size0 = Pre_Subreg_info[1];

        for (int r = 0; r < Pre_Num_Subreg && Pre_Subreg_info[r * record_size0 + 2] == 1; ++r) {
            Pre_FinG_Num++;
        }
        Pre_EMDF_Num = Pre_Num_Subreg - Pre_FinG_Num;
    }

    int FinG_Num = 0;
    int EMDF_Num = 0;
    int record_size = 0;

    if (SubregInfo_in != nullptr) {
        int Num_Subreg = SubregInfo_in[0];
        record_size = SubregInfo_in[1];

        for (int r = 0; r < Num_Subreg && SubregInfo_in[r * record_size + 2] == 1; ++r) {
            FinG_Num++;
        }
        EMDF_Num = Num_Subreg - FinG_Num;
    }

    // Poro volume
    double* Pre_Porv = new double[Pre_Net_size];
    if (a_Pre_Porv.length == Pre_Net_size) {
        for (int i = 0; i < Pre_Net_size; ++i)
            Pre_Porv[i] = a_Pre_Porv.fdata[i];
    }
    else {
        for (int i = 0; i < Pre_Net_size; ++i) 
            Pre_Porv[i] = 1.0;
    }

    int* Pre_G2N_map = new int[Pre_Gross_size];
    for (int i = 0; i < Pre_Gross_size; ++i) {
        Pre_G2N_map[i] = origin_Gross2Net[i];
    }

    /************************************************************************/
    /*          Сд��ѭ������������Ч������, ��д��ѭ�����������ܱ��        */
    /************************************************************************/

    /* Step 1: Host_grid��"ԭ"���������ӳ�� */
    UGSIM::sp_vec<vector<int>> Host2GrossFinG(Base_NLv0* PoroLvl);
    UGSIM::sp_vec<vector<int_doub>> Host2NetFinG(Base_NLv0* PoroLvl);

    if (Pre_Gross_fin_size > 0 && Pre_Host_map != nullptr) {

        for (int lvl = 0; lvl < PoroLvl; ++lvl) {
            int START = Base_NLv0 + Pre_Total_NLv0 * lvl;
            int END = Pre_Total_NLv0 * (lvl + 1);

            for (int I = START; I < END; ++I) {
                int HOST = Pre_Host_map[I];
                Host2GrossFinG[HOST].push_back(I);

                int i = Pre_G2N_map[I];
                if (i != iDead) {
                    Host2NetFinG[HOST].push_back(int_doub(i, Pre_Porv[i]));
                }
            }
        }

        int Host_size = Base_NLv0 * PoroLvl;
        for (int I = 0; I < Host_size; ++I) {
            if (Host2NetFinG.is_filled(I)) {
                vector<int_doub>& Weights = Host2NetFinG[I];
                double Porv = 0.0;

                for (int j = 0; j < Weights.size(); ++j) {
                    Porv += Weights[j].db;
                }
                for (int j = 0; j < Weights.size(); ++j) {
                    Weights[j].db /= max(Porv, InfSmall);
                }
            }
        }
    }

    /* �������� */
    Moves.reserve(Net_size);

    /* Step 2: ��������"ԭ"����ӳ��, push��Moves */
    int FLvlNum = FLvlBaseNum + FLvlFinNum;
    int MLvlNum = MLvlBaseNum + MLvlFinNum;
    for (int lvl = 0; lvl < PoroLvl; ++lvl) {
        int start = (lvl == 0) ? 0 : FLvlNum + (lvl - 1) * MLvlNum;
        int end = (lvl == 0) ? FLvlBaseNum : FLvlNum + (lvl - 1) * MLvlNum + MLvlBaseNum;
        int HOST_TO_PREG = Pre_Gross_fin_size / PoroLvl * lvl;

        for (int to = start; to < end; ++to) {
            int HOST = (Host_map_in != nullptr) ? Host_map_in[Net2Gross_in[to]] : Net2Gross_in[to];
            int from = Pre_G2N_map[HOST + HOST_TO_PREG];

            if (from != iDead) {
                Moves.push_back(transfer(from, to, 1.0));
            }
            else {
                if (Host2NetFinG.is_filled(HOST)) {
                    const vector<int_doub>& Weights = Host2NetFinG[HOST];
                    for (int j = 0; j < Weights.size(); ++j) {
                        Moves.push_back(transfer(Weights[j].i, to, Weights[j].db));
                    }
                }
                else {
                    Moves.push_back(transfer(-1, to, 1.0));
                }
            }
        }
    }

    /* Step 3: ��������--"ԭ"��������ӳ�� */
    map<string, int> Name_to_Carfin;
    vector<int> Pre_FinG_START;

    if (Pre_FinG_Num > 0 && Pre_Subreg_info != nullptr) {
        int FinG_START = Base_NLv0;

        for (int r = 0; r < Pre_FinG_Num; ++r) {
            const int* Pre_Subreg_data = Pre_Subreg_info + r * record_size0 + 2;
            char ctmp[MAX_STRING_LEN] = "";
            ia_to_str(Pre_Subreg_data + REG_INFO_SIZE, MAX_STRING_LEN, ctmp);
            Name_to_Carfin[string(ctmp)] = r;

            Pre_FinG_START.push_back(FinG_START);
            FinG_START += Pre_Subreg_data[1] * Pre_Subreg_data[2] * Pre_Subreg_data[3]; // NX*NY*NZ
        }
        Pre_FinG_START.push_back(FinG_START);
    }

    map<int, int> FinG_to_FinG;

    if (FinG_Num > 0 && SubregInfo_in != nullptr && Host_map_in != nullptr) {

        for (int lvl = 0; lvl < PoroLvl; ++lvl) {
            int PREG_OFFSET = lvl * Pre_Total_NLv0;
            int HOST_TO_PREG = Pre_Gross_fin_size / PoroLvl * lvl;
            int TO = Base_NLv0 + lvl * Total_NLv0;

            for (int r = 0; r < FinG_Num; ++r) {
                const int* Subreg_data = SubregInfo_in + r * record_size + 2;
                char ctmp[MAX_STRING_LEN] = "";
                ia_to_str(Subreg_data + REG_INFO_SIZE, MAX_STRING_LEN, ctmp);

                int r0 = 0;
                bool found = Key_to_Ind(Name_to_Carfin, string(ctmp), r0);

                bool match = false;
                if (found && Pre_Subreg_info != nullptr) {
                    const int* Pre_Subreg_data = Pre_Subreg_info + r0 * record_size0 + 2;
                    match = (ia_compare(Subreg_data, Pre_Subreg_data, REG_INFO_SIZE) == 0);
                }

                int BOX_SIZE = Subreg_data[1] * Subreg_data[2] * Subreg_data[3];

                if (match) { // ���һ: ��ȫƥ��, һһӳ��
                    int START = Pre_FinG_START[r0] + PREG_OFFSET;

                    for (int I = 0; I < BOX_SIZE; ++I) {
                        int to = Gross2Net_in[TO];

                        if (to != iDead) {
                            int from = Pre_G2N_map[I + START];
                            if (from != iDead) {
                                Moves.push_back(transfer(from, to, 1.0));
                            }
                            else { // mismatch
                                Moves.push_back(transfer(-1, to, 1.0));
                            }
                        }

                        if (lvl == 0) {
                            FinG_to_FinG[TO] = I + START;
                        }

                        ++TO;
                    }
                }
                else { // �����: ����ȫƥ��, ��ԭ��������ȡֵ               
                    for (int I = 0; I < BOX_SIZE; ++I) {
                        int to = Gross2Net_in[TO];

                        if (to != iDead) {
                            int HOST = Host_map_in[TO];
                            int from = Pre_G2N_map[HOST + HOST_TO_PREG];

                            if (from != iDead) {
                                Moves.push_back(transfer(from, to, 1.0));
                            }
                            else {
                                if (Host2NetFinG.is_filled(HOST)) {
                                    const vector<int_doub>& Weights = Host2NetFinG[HOST];
                                    for (int j = 0; j < Weights.size(); ++j) {
                                        Moves.push_back(transfer(Weights[j].i, to, Weights[j].db));
                                    }
                                }
                                else { // mismatch
                                    Moves.push_back(transfer(-1, to, 1.0));
                                }
                            }
                        }

                        ++TO;
                    } // for I 
                } // if not match
            } // for r
        } // for lvl
    }

    // Step 4: EMDF--"ԭ"EMDFӳ��
    map<string, int> Name_to_EMDF;
    vector<int> Pre_Emdf_START;

    if (Pre_Subreg_info != nullptr) {
        int EMDF_START = Pre_Gross_size - Pre_Gross_emdf_size;

        for (int r = Pre_FinG_Num; r < Pre_FinG_Num + Pre_EMDF_Num; ++r) {
            const int* Pre_Subreg_data = Pre_Subreg_info + r * record_size0 + 2;
            char ctmp[MAX_STRING_LEN] = "";
            ia_to_str(Pre_Subreg_data + REG_INFO_SIZE, MAX_STRING_LEN, ctmp);
            Name_to_EMDF[string(ctmp)] = r - Pre_FinG_Num;

            Pre_Emdf_START.push_back(EMDF_START);
            EMDF_START += Pre_Subreg_data[1];  // face num
        }
        Pre_Emdf_START.push_back(EMDF_START);
    }

    if (EMDF_Num > 0 && SubregInfo_in != nullptr && Host_map_in != nullptr) {

        int PREG_OFFSET = (PoroLvl == 1) ? 0 : Pre_Total_NLv0; // ʹ"ԭ"��������ƫ�Ƶ�Lvl1����
        int HOST_OFFSET = (PoroLvl == 1) ? 0 : Base_NLv0;      // ʹHost2NetFinGƫ�Ƶ�Lvl1����
        int TO = Gross_size - Gross_emdf_size;                 // TO���ѷ��������, �ܱ��

        for (int r = FinG_Num; r < FinG_Num + EMDF_Num; ++r) {
            const int* Subreg_data = SubregInfo_in + r * record_size + 2;
            char ctmp[MAX_STRING_LEN] = "";
            ia_to_str(Subreg_data + REG_INFO_SIZE, MAX_STRING_LEN, ctmp);

            int r0 = 0;
            bool name_match = Key_to_Ind(Name_to_EMDF, string(ctmp), r0);

            map<int, int> Pre_BG_emdf;  // BgGrid to fracture grid map
            if (name_match) {
                int START = Pre_Emdf_START[r0];
                int END = Pre_Emdf_START[r0 + 1];

                for (int I = START; I < END; ++I) {
                    int i = Pre_G2N_map[I];
                    int BG = Pre_Host_map[I];
                    if (i != iDead) {
                        Pre_BG_emdf[BG] = i;
                    }
                }
            }

            int FACE_SIZE = Subreg_data[1];
            for (int I = 0; I < FACE_SIZE; ++I) {
                int to = Gross2Net_in[TO];

                if (to != iDead) {
                    int BG = Host_map_in[TO]; // back ground grid

                    if (BG < Base_NLv0) { // δ����EMDF
                        bool match = false;

                        if (name_match) {                    // ���ѷ�����ȡֵ
                            int face = 0;
                            if (Key_to_Ind(Pre_BG_emdf, BG, face)) {  // 1. δ����EMDF --> δ����EMDF 
                                Moves.push_back(transfer(face, to, 1.0));
                                match = true;
                            }
                            else {                           // 2. ����EMDF --> δ����EMDF 
                                if (Host2GrossFinG.is_filled(BG)) {
                                    const vector<int>& FinBGs = Host2GrossFinG[BG];
                                    vector<int_doub> Weights;

                                    for (int j = 0; j < FinBGs.size(); ++j) { // �ռ�sub_face
                                        int sub_face = 0;
                                        if (Key_to_Ind(Pre_BG_emdf, FinBGs[j], sub_face)) {
                                            Weights.push_back(int_doub(sub_face, Pre_Porv[sub_face]));
                                        }
                                    }

                                    if (!Weights.empty()) {
                                        double Porv = 0.0;
                                        for (int j = 0; j < Weights.size(); ++j) {
                                            Porv += Weights[j].db;
                                        }
                                        for (int j = 0; j < Weights.size(); ++j) {
                                            Moves.push_back(transfer(Weights[j].i, to, Weights[j].db / Porv));
                                        }
                                        match = true;
                                    }
                                }
                            }
                        }

                        if (!match) {                       // �ӻ�������ȡֵ
                            int from = Pre_G2N_map[BG + PREG_OFFSET];
                            if (from != iDead) {            // 3.δ���ܻ������� --> δ����EMDF
                                Moves.push_back(transfer(from, to, 1.0));
                                match = true;
                            }
                            else {                          // 4.���ܻ������� --> δ����EMDF
                                if (Host2NetFinG.is_filled(BG + HOST_OFFSET)) {
                                    const vector<int_doub>& Weights = Host2NetFinG[BG + HOST_OFFSET];
                                    for (int j = 0; j < Weights.size(); ++j) {
                                        Moves.push_back(transfer(Weights[j].i, to, Weights[j].db));
                                    }
                                    match = true;
                                }
                            }
                        }

                        if (!match) { // mismatch
                            Moves.push_back(transfer(-1, to, 1.0));
                        }
                    }
                    else {  // ����EMDF
                        bool match = false;
                        int FinG0 = 0;
                        bool Found_FinG = Key_to_Ind(FinG_to_FinG, BG, FinG0);

                        if (name_match) {                   // ���ѷ�����ȡֵ
                            if (Found_FinG) {               // 1.����EMDF --> ����EMDF, һ��һ
                                int face = 0;
                                if (Key_to_Ind(Pre_BG_emdf, FinG0, face)) {
                                    Moves.push_back(transfer(face, to, 1.0));
                                    match = true;
                                }
                            }
                            else {
                                int BASE = Host_map_in[BG];
                                int face = 0;
                                if (Key_to_Ind(Pre_BG_emdf, BASE, face)) {  // 2.δ����EMDF --> ����EMDF
                                    Moves.push_back(transfer(face, to, 1.0));
                                    match = true;
                                }
                                else {                       // 3.����EMDF --> ����EMDF, ���һ
                                    if (Host2GrossFinG.is_filled(BASE)) {
                                        const vector<int>& FinBGs = Host2GrossFinG[BASE];
                                        vector<int_doub> Weights;

                                        for (int j = 0; j < FinBGs.size(); ++j) { // �ռ�sub_face
                                            int sub_face = 0;
                                            if (Key_to_Ind(Pre_BG_emdf, FinBGs[j], sub_face)) {
                                                Weights.push_back(int_doub(sub_face, Pre_Porv[sub_face]));
                                            }
                                        }

                                        if (!Weights.empty()) {
                                            double Porv = 0.0;
                                            for (int j = 0; j < Weights.size(); ++j) {
                                                Porv += Weights[j].db;
                                            }
                                            for (int j = 0; j < Weights.size(); ++j) {
                                                Moves.push_back(transfer(Weights[j].i, to, Weights[j].db / Porv));
                                            }
                                            match = true;
                                        }
                                    }
                                }
                            }
                        }

                        if (!match) {                       // �ӻ�������ȡֵ                   
                            if (Found_FinG) {               // 4.���ܻ������� --> ����EMDF, һ��һ
                                int from = Pre_G2N_map[FinG0 + PREG_OFFSET];
                                if (from != iDead) {
                                    Moves.push_back(transfer(from, to, 1.0));
                                    match = true;
                                }
                            }
                            else {
                                int BASE = Host_map_in[BG];
                                int from = Pre_G2N_map[BASE + PREG_OFFSET];

                                if (from != iDead) {        // 5.δ���ܻ������� --> ����EMDF
                                    Moves.push_back(transfer(from, to, 1.0));
                                    match = true;
                                }
                                else {                      // 6.���ܻ������� --> ����EMDF, ���һ
                                    if (Host2NetFinG.is_filled(BASE + HOST_OFFSET)) {
                                        const vector<int_doub>& Weights = Host2NetFinG[BASE + HOST_OFFSET];
                                        for (int j = 0; j < Weights.size(); ++j) {
                                            Moves.push_back(transfer(Weights[j].i, to, Weights[j].db));
                                        }
                                        match = true;
                                    }
                                }
                            }
                        }

                        if (!match) {
                            Moves.push_back(transfer(-1, to, 1.0));
                        }
                    }
                }

                ++TO;
            }

        }
    }

    delete[] Pre_Porv;
    delete[] Pre_G2N_map;

    // Step 5: ����
    int m0 = 0;
    while (m0 < Moves.size() && Moves[m0].from == -1) { ++m0; };

    if (m0 < Moves.size()) {
        // ������в�ƥ�䣬������������Ԫ���ֲ� 
        int mis_match = m0;

        for (int m = 0; m < m0; ++m) {
            Moves[m].from = Moves[m0].from;
            Moves[m].mult = 1.0;
        }

        for (int m = 1; m < Moves.size(); ++m) {
            if (Moves[m].from == -1) {
                Moves[m].from = Moves[m - 1].from;
                Moves[m].mult = 1.0;
                ++mis_match;
            }
        }

#if _DEBUG
        // Check unity
        {
            double* work = new double[Net_size];
            for (int i = 0; i < Net_size; ++i) {
                work[i] = 0.0;
            }

            int exceeds = 0;
            for (vector<transfer>::const_iterator it = Moves.begin(); it != Moves.end(); ++it) {
                if (it->to < Net_size && it->from < Pre_Net_size) {
                    work[it->to] += it->mult;
                }
                else {
                    ++exceeds;
                }
            }

            int non_unit = 0;
            for (int i = 0; i < Net_size; ++i) {
                if (fabs(work[i] - 1.0) > NearZero) {
                    ++non_unit;
                }
            }
            delete[] work;

            if (exceeds > 0) {
                cout << " Error: transfer sequence exceeds grid index range" << endl;
            }
            if (non_unit > 0) {
                cout << " Error: transfer sequence does not exactly sum up to unit" << endl;
            }
        }
#endif

        // Step 6: �滻����ṹ
        for (int i = 0; i < 32; ++i) {
            GridInfo[i] = GridInfo_in[i];
        }

        Gross2Net.resize(Gross_size);
        for (int i = 0; i < Gross_size; ++i) {
            Gross2Net[i] = Gross2Net_in[i];
        }

        Net2Gross.resize(Net_size);
        for (int i = 0; i < Net_size; ++i) {
            Net2Gross[i] = Net2Gross_in[i];
        }

        return mis_match;
    }
    else {
        Moves.clear();
        return -2;
    }
};